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Abstract 


Signal-processing on graphs has developed into a very active field of research 
during the last decade. In particular, the number of applications using frames con¬ 
structed from graphs, like wavelets on graphs, has substantially increased. To at¬ 
tain scalability for large graphs, fast graph-signal filtering techniques are needed. 

In this contribution, we propose an accelerated algorithm based on the Lanczos 
method that adapts to the Laplacian spectrum without explicitly computing it. The 
result is an accurate, robust, scalable and efficient algorithm. Compared to existing 
methods based on Chebyshev polynomials, our solution achieves higher accuracy 
without increasing the overall complexity significandy. Furthermore, it is particu¬ 
larly well suited for graphs with large spectral gaps. □ 

Index terms — Spectral graph theory, graph signal-processing, graph filter, Chebyshev 
polynomial, Lanczos method 

1 Introduction 

Graphs conveniently represent data on irregular geometric structures as they arise in 
numerous application domains such as social, energy, transportation or neuronal net¬ 
works. In all of these fields, different pieces of information are connected to each other 
and these connections can be modelled by a graph. For every link, an edge is drawn 
with an associated weight that represents the similarity between the two elements (ver¬ 
tices) it connects. For example, in a sensor network acquiring environmental mea¬ 
surements, it makes sense to choose the weight inversely proportional to the physical 
distance. A data signal lives on nodes and can be visualized as a collection of samples. 
We refer to those as graph-signals, see Figure[T]for an example. 

This framework, called graph signal-processing, has recently become a general 
field of research ifT^ flTll and applications of graph signal-processing can be found in 
many different areas. In machine vision, automatic text classification or more generally 
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Figure 1: Synthesized sensor network. Colored circles represent vertices and gray links 
represent edges. Values of the signal are displayed with different colors. 

machine learning, graphs are used to represent similarities between data points and re¬ 
sult in algorithms for semi-supervised learning ll24l[T9l [n i23]l . Graph signal-processing 
is often used to leverage intrinsic links when processing data. For instance, many ap¬ 
plications in image processing benefit from the use of graphs that describe connections 
between local patches in the image (e.g. lfT3ll22l and references therein). 

Spectral graph filtering is one of the basic building blocks in the applications dis¬ 
cussed above. Classical filtering techniques rely on the Fourier transform, which can 
be done inexpensively with a cost of 0{N log(iV)). Extending this transform to graphs 
requires the diagonalisation of the graph Laplacian, which becomes prohibitively ex¬ 
pensive for larger graphs. To overcome this issue, an efficient approximate method to 
perform spectral graph filtering was proposed in ifTOll . This method, based on Cheby- 
shev polynomial approximation, only requires multiple application of the Laplacian 
operator and thus leads to a scalable, efficient memory-saving and error-controlled al¬ 
gorithm. 

Instead of Chebyshev polynomial approximation, we propose to use the Lanczos 
method to perform spectral graph filtering. The resulting algorithm is also scalable, in 
addition to being superior in accuracy to its predecessor. 

The rest of this article is organized as follows. In Section 2, we summarize the 
basics of signal-processing on graphs, including the definition of graph filtering. We 
also recall the method of Hammond et al. Qo). Section 3 describes the Lanczos method 
and its application to graph filtering. Numerical experiments are presented in Section 
4. 
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2 Signal-processing on graphs 

2.1 Graph nomenclature 

We consider a weighted undirected graph Q = {V, £, W} with a set of vertices V, a 
set of edges £, and a weight function W ; V x V —R. The vertices are indexed from 
1,..., = I V| and each entry of the weight matrix W G contains the weight 

of the edge connecting the corresponding vertices; Wij = yV{vi,Vj). If there is no 
edge between two vertices, the weight is set to 0. It is assumed that Wij = Wj^i, that 
is, is a symmetric matrix. For a vertex Vi € V, the degree d{i) is defined as the sum 
of the weights of incident edges: d{i) = 

In this framework, a graph signal is defined as a function s : V —R assigning a 
value to each vertex. It is convenient to consider a signal s as a vector of size N with 
the i*** component representing the signal value at the i*** vertex. 

One of the most fundamental concepts for weighted undirected graphs is the (com¬ 
binatorial) graph Laplacian £ defined as C = D — W, where D is the diagonal degree 
matrix with diagonal entries Du = d{i). Alternative definitions of graph Laplacians 
include the normalized Laplacian £„ = D^LD^ = {D — W)D'^. For simplicity, 
we restrict ourselves to £, but the method presented in this paper easily extends to other 
Laplacians. 

The Laplacian C is always symmetric positive semi-definite and can thus be de¬ 
composed as 

C = UKU\ 

where U = [uq, ... ,un-i] is an orthogonal matrix and U* denotes its transpose. 
Without loss of generality, we order the set of eigenvalues as follows: 0 = Aq < Ai < 
... < \ n -1 = Amaxi see im for more details on spectral graph theory. The matrix U 
defines the graph Fourier basis Efiiini, leading to the graph Fourier transform s = U*s 
and its inverse s = Us. 

2.2 Graph filters 

In the classic setting, applying a filter to a signal is performed by convolution, which 
corresponds to point-wise multiplication in the spectral domain. Similarly, filtering a 
graph-signal is peformed by multiplication with a filter in the graph Fourier domain. 
A graph filter is defined by a continuous function g : R+ R. To obtain its discrete 
coefficients, this function is evaluated at each eigenvalue: g{Xe) for £ = 0,..., N — 1. 
The filtering operation then corresponds to s'{i) = g{Xe) ■ s{£), where s' is the filtered 
signal. Equivalently, using matrix notation, we have 

s' = Ug{A)U*s, (1) 

where g{A) is the diagonal matrix containing the coefficients g{Xg) on the diagonal. 
In terms of matrix functions ini, the relation ([T]i can be compactly expressed as s' = 
g{C)s with g{C) := Ug{A)U*. 
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2.3 Fast filtering via Chebyshev polynomials 

The graph filtering operation described above is based on the graph Fourier transform. 
Unfortunately, the graph Fourier basis needed for performing this transform requires 
the diagonalization of the graph Laplacian, which takes 0{N^) operations and 0{N'^) 
memory when using standard techniques. This is feasible for graphs with only a few 
thousand vertices. To be able to tackle problems of larger size, more efficient methods 
are needed, one of which is presented in (nni and summarized in this section. 

Filtering in the vertex domain. To avoid the Fourier transform, we perform the 
filtering operation in the vertex domain using only the Laplacian operator. Applying 
this operator corresponds to multiplying the signal in the spectral domain with the 
eigenvalues: 


Cs = As. 


This is equivalent to filtering with g{x) = x. Using this relation recursively and ex¬ 
ploiting linearity, we can apply any polynomial filter g(x) = oq + aix + • • • + qkx^ 
to a signal s with the following formula: 


s' = Ug{K)U*s = {aol + aiC + --- + s . 


( 2 ) 


Chebyshev polynomial approximation. The ability to apply polynomial filters 
efficiently suggests to approximate a given filter function with a suitable polynomial. 
For approximating functions on real intervals, Chebyshev polynomials are usually the 
preferred choice because of numerical stability considerations and the fact that they can 
be evaluated efficiently by three-term recurrences. We refer to mni for a more detailed 
discussion on the choice of Chebyshev polynomials in signal-processing on graphs and 
to, e.g., mi for an introduction to polynomial approximation. 

The mth Chebyshev polynomial Tm{y) is generated using the recurrence relation 
Tm{y) = 2xTm-i{y) - Tm-i{y) with To{y) = 1 and Ti(y) = y. For y £ [-1,1], 
these polynomials possess the following well-known properties: 

1. they admit the closed form expression Tm{y) = cos(m arccos(t/)); 

2. they are bounded, i.e., Tm(y) £ [—1,1]; 




admits a 


convergent Chebyshev series 



m—1 


with the Chebyshev coefficients 
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Since our filter g is evaluated on the eigenvalues of the Laplacian, we need to map the 
interval [—1,1] to the interval [0, Amax] using the transformation x = {y + 1). 

Defining Tm{x) = Tm -we obtain 


9[x) 


1 

2 


OO 

Co ^ ^ ^ 

m—1 


(3) 


forx e [0, Amax], with 

Cm = — j cos{m6)g (cos(0) + 1)^ d9. 

Fast filtering algorithm. At this point, we can derive the iterative algorithm for 
filtering a signal s with g. The recurrence relation for the transformed Chebyshev 

polynomials becomes fm{x) = 2 - Tm-i{x) — Tm- 2 {x). On the matrix 

level, this yields, using (|2]i; 

fm{C)s = 2 f ^ - /) fm-l{C)s - fm-2{C)s. 

\ ^max / 

Combined with this finally leads to the following expression for filtering a signal 
5 : 

^ 00 

s' = g{C) = -^cqIs + CmTm{C)s. 

m—1 

When implemented, we truncate this sum at a defined order M. Assuming that \S\ > 
N, the computational cost of this algorithm scales linearly with the number of edges 
0{M\£\). In most applications, the Laplacian is sparse, \£\ ^ which results in 
a fast algorithm. Moreover, apart from storing the Laplacian, the additional memory 
consumed by this algorithm is only 4A^. 


3 Accelerated filtering using Lanczos 

Given the graph Laplacian C G and a nonzero vector s € R^, the Lanc¬ 

zos method 13 shown in Algorithm [T] below computes an orthonormal basis Vm = 
[ui,..., Vm] of the Krylov subspace Km{£, s) = span{s, £s,..., The 

computational cost of Algorithm is 0{M ■ \£\). The storage of the basis Vm requires 
MN additional memory, which can be avoided using two passes of the algorithm or 
restart techniques; see, e.g., Q for more details. 

The scalars produced by Algorithm[T]can be arranged into a symmetric tridiagonal 


matrix Hm G satisfying 

ai 

P 2 




P 2 

0-2 



V^iCVm =Hm = 


/33 

as 

Pm 





Pm olm 
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Algorithm 1 Lanczos method 

Input: Symmetric matrix C € vector s 0, M G N. 


Output: Vm = [ft) ■ • ■, vm] with orthonormal columns, scalars ai,. 

., aM C R and 


P 2 , ■ ■ ■, Pm S M. 


1 

Vl ^ s/||s ||2 


2 

for j = 1,2,..., M do 


3 

S 

II 

bi 


4 

Oij = v*w 


5 

Vj + l = w- VjUj 


6 

if j > 1 then 


7 

Vj -(_ 1 i '^j — 1 l^j — 1 


8 

end if 


9 

Pj = l%-i-i I 2 


10 

Vj + i = Vj + i/Pj 


11 

end for 



In floating-point arithmetics, the orthogonality of the basis produced by Algorithm [T] 
may get quickly lost and reorthogonalization is needed El. 

Given a continuous function g : [0,Amax] —R and a vector s, the following 
approximation to p(£)s was proposed by Gallopoulos and Saad in 0: 

g{C)s \\s\\ 2 VMgiHM)ei := gn, (4) 

where ei € is the first unit vector. Because of eigenvalue interlacing, the eigen¬ 
values of Hm are contained in the interval [0, Amax] and hence the expression g{HM) 
is well-defined. Typically, M N, rendering the evaluation of g{HM) inexpensive. 
The overall cost of our Lanczos-based approximation of graph-signal filtering, which 
consists of applying Algorithm [T] and evaluating (|4|i, is therefore between 0{M ■ |) 

and 0{M ■ \S\ -f M'^N), depending on how the reorthogonalization is performed. 

Approximation quality. The following theorem provides some insight into the ac¬ 
curacy of the approximation gM obtained by the Lanczos method. 

Theorem 1 (||9] Corollary 3.4]). Let C € be symmetric with eigenvalues con¬ 

tained in the interval [0, Amax] ond let g : [0, Amax] —>■ M be continuous. Then 

||p(£)s-pm||2 < 2||s|| 2- min max \g{z) - p{z)\, 

P&Vm^i zG[0,A„ax] 

where Vm-i denotes all polynomials of degree at most M — 1. 

Theorem [T] shows that the error is bounded by the best polynomial approxima¬ 
tion IH of g on [OjAmax]. Compared to the Chebyshev approximation from Sec¬ 
tion JZS] this shows that - up to a multiple of two - the Lanczos-based approximation 
gM can be expected to provide at least the same accuracy. 

However, the Lanczos-based approximation can sometimes be expected to perform 
significantly better because of its ability to adapt to the eigenvalues of C. This phe¬ 
nomenon is well-understood for Krylov subspace approximations to solutions of linear 
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systems (see, e.g., (H Section 3.1]) and extends to matrix functions as well. To illus¬ 
trate such a situation, suppose that the eigenvalue 0 of £ is well separated from the 
other eigenvalues contained in the (smaller) interval [Ai, Amax]- Then the eigenvalues 
of Hm can be expected to exhibit a similar behavior IfTSi Chapter 6]. 

By choosing a polynomial that interpolates the eigenvalue 0 exactly and approxi¬ 
mates g on [Ai, Amax], the bound of Theorem[T]can approximately be replaced by 

||p(£)s-pm||2 < 2||s||2- min max \g{z) - p{z)\. 

P^PM-2 zelAijAmax] 

At the expense of losing a degree of freedom in the choice of the polynomial, the width 
of the approximation interval becomes smaller, which in turn leads to significantly 
improved approximation rates lfT4ll . In contrast, the Chebyshev approximation cannot 
adapt to the eigenvalues of £ and hence its approximation rate is driven by the larger 
interval [0, Amax]- 

Stopping criteria. Ideally, M should be chosen to be the smallest integer such that 
IjeMlb = 1|5(£)s-5m 1|2 < e 

holds for some prescribed accuracy e > 0. As proposed in ll20l . we estimate cm as the 
difference of two (consecutive) approximations: 

e-M ~ 9M+i — givi (5) 


for a small value of ). 

4 Numerical experiments 

In this section we present numerical experiments to demonstrate the numerical be¬ 
haviour of the Lanczos method and compare it with the Chebyshev method. All ex¬ 
periments were performed with the signal processing toolbox GSPBox IfT^ and can be 
downloaded at https : / /lts 2 . epf 1 . ch/lanczos-f iltering/ 

Example 1. We first consider the sensor graph and the Erdos-Renyi random graph ID, 
with N = 500 nodes. In the latter case, each edge is included in the graph with 
probability p = 0.04. As a filterbank, we use a collection of translated windows 
g{t) = sin(0.57r cos(7rf)^)x[_i/2.i/2], where xi is the characteristic function of 1 1211 . 
We choose to adapt the filterbank to the graph spectrum, since non-adapted filters lead 
to very coherent frames for graphs with eigenvalues not uniformly spread along the 
spectrum US). 

First, we test the reliability of the error estimate (H) for j = 3. Figure |2] confirms 
that the estimate UffM-i-s ~ <?m 1|2 closely follows the true error ||eMl|2- 

In Figures[3]and|5]we observe that the Fanczos method yields better approximation, 
especially in the case when the spectrum of graph Faplacian is not uniformly distributed 
and has a large (relative) spectral gap, as predicted in Section[3 

Example 2. The purpose of this example is to investigate how the approximation 
by the Chebyshev polynomial method and the Fanczos method depends on the edge 
sparsity of the Erdos - Renyi graph. We consider graphs with N = 1000 nodes. 
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Figure 2; Sensor graph (left) and Erdos - Renyi graph (right). Comparison of the 
approximation error in the Lanczos method and the estimate Q. 



Figure 3: Sensor graph (left) and Erdos - Renyi graph (right). Comparison of the 
approximation error using the Chebyshev method and the Lanczos method with respect 
to the order of approximation. 


The order of Chebyshev polynomial and Lanczos method are set to M = 30. As a 
filterbank, we the use mexican hat wavelet, with a mother window 


gh[>^t) = Xi ■ exp(-A|), 

with Xg eigenvalues of graph Laplacian. In order to obtain a complete filterbank, we 
add a low pass filter gi{\t) = exp(—A^). 

We notice that the Lanczos method exhibits excellent numerical behaviour in com¬ 
parison with the Chebyshev method, particularly when using non-adapted filters. As 
the probability p increases, the (relative) spectral gap increases as well, leading to more 
accurate Lanczos approximation (see Figure|6]l. 

5 Conclusion 

In this letter, a scalable, efficient and accurate method to perform signal filtering on 
graphs based on the Lanczos method is proposed. We compare it to the existing method 
based on the Chebyshev polynomial approximation. The advantage of filtering using 
the Lanczos approximation is that no a priori information about the spectrum of C is 
needed. In all numerical experiments the proposed method outperforms the Chebyshev 
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Figure 4: Sensor graph (left) and Erdos - Renyi graph (right). Filterbanks adapted to 
spectra of graph Laplacians. 
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Figure 5: Sensor graph (left) and Erdos - Renyi graph (right). Distribution of graph 
Faplacians’ eigenvalues. 




Figure 6: Comparison of the approximation error using the Chebyshev method and 
the Fanczos method with respect to p, with non-adapted filterbank (left) and adapted 
filterbank (right). 


approach. That is especially visible in case of a non-adapted filterbank. Futhermore, 
the Fanczos approximation is superior to the Chebyshev approximation when filtering 
on random graphs with well separated eigenvalues. 
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Figure 7: Non-adapted filterbank (left) and adapted filterbank (right). 
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